Self-Organized Criticality and Thermodynamic formalism. 
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We introduce a dissipative version of the Zhang's model of Self-Organized Criticality, where a 
parameter allows to tune the local energy dissipation. We analyze the main dynamical features 
of the model and relate in particular the Lyapunov spectrum with the transport properties in the 
stationary regime. We develop a thermodynamic formalism where we define formal Gibbs measure, 
partition function and pressure characterizing the avalanche distributions. We discuss the infinite 
size limit in this setting. We show in particular that a Lee- Yang phenomenon occurs in this model, 
for the only conservative case. This suggests new connexions to classical critical phenomena. 
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c/^ . 

I In 1988, Bak, Tang and Wiesenfeld (BTW) [1] proposed for the first time a mechanism allowing a dynamical 
system to reach "spontaneously" a steady state exhibiting some features of a thermodynamic system at a critical 
^— ^ • point. Namely, by its only internal reorganization in reaction to (stationary) external world perturbations, the system 
• reaches a stationary state with power law statistics. This effect, called by these authors Self-Organized Criticality 
^ ' (SOC), was quite unexpected since attaining the critical state of a thermodynamic system usually needs a fine tuning 
(-H , of some control parameter (temperature, magnetic field, etc..) that is at first sight absent from the definition of the 
BTW model and its many variants [2,3]. This is a first difference between the SOC systems and thermodynamic 
^ ' systems exhibiting a second order phase transition. Furthermore, at stationarity, there is a constant flux through 
I the system since the stationary flux of external perturbations is dissipated in the bulk or at the boundaries. This 
, corresponds to a non-equilibrium situation and it is therefore believed that usual equilibrium statistical mechanics 
r treatments using the concept of Gibbs measure, free energy, etc ... cannot be applied in this case. Moreover, there 

(N 

is an implicit notion of "thermodynamic limit" in SOC systems in order to reach the critical state. However, the 
^ ' absence of a thermodynamic setting raises difficulties in establishing a proper definition of the thermodynamic limit 
(-H ] and only a few paper (like [4]) have focused on this aspect. Consequently, the results concerning this limit rely mainly 
^ : on extrapolations from the finite size numerical data and the properties of the limiting system have to be conjectured 
from the finite size dynamical system. Unfortunately, these numerical results must be handled properly. This implies 
in particular to have the right finite size scaling form and a fair control of the bias induced by numerics. 

In [5-8] we proposed an analysis of one SOC model, the Zhang's model [9], using the tools and concepts from 
dynamical system theory and ergodic theory. In particular, our results opened the perspective to construct a thermo- 
dynamic formalism in the sense of Sinai, Ruelle, Bowen [10-12]. In this setting, one defines in particular an extension 
of the notion of Gibbs measure, where the Hamiltonian is replaced by a dynamically relevant quantity. We discussed 
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in particular in [5,7] the possibility to use this formalism to relate the microscopic dynamics to the macroscopic 
characteristics of the critical state. Actually, non trivial results relating the fractal structure of the attractor and the 
probability distribution of avalanche sizes were obtained. Moreover, we argued that this formalism could be used to 
define properly the stationary state of the finite size dynamical system, and then an extrapolation of its properties 
when the size tends to infinity could be done. 

This paper and [13] are intending to develop this aspect. More precisely, our aims are: 

• Construct the equivalent of finite volume Gibbs measures (in the sense of Sinai, Ruelle, Bowen) in order to 
have the equivalent of the statistical mechanics generating functions : partition function and free energy, to 
characterize the stationary state. Since there is no Hamiltonian in the dynamics this construction must be done 
in a more general setting. 

• Link these generating functions to dynamical properties. 

• Relate these functions to the avalanche distributions. 

• Analyze the behavior of these functions for finite size and attempt to characterize the SOC state by extrapolating 
to the limit. 

• Have a theoretical control of the validity of the extrapolations made from numerical simulations. 

The two last points are the main topic of a separated and independent paper [13]. In the present paper we will 
mainly focus on the three first aspects. For, we introduce a non conservative version of the Zhang model where 
a parameter h controls the local energy dissipation. This parameter acts somehow as a temperature and we show 
that it has to be tuned to ft, = in order to obtain a critical state in the thermodynamic limit. Consequently, our 
model "self-organizes" into a "critical" state only in the conservative case. However, the presence of this additional 
parameter gives us some additional freedom to control the dynamics, and allows us to define an additional critical 
exponent. Moreover, it allows us to outline the role of boundary dissipation in the process of self-organization into a 
critical state. 

The paper is structured as follows. In the first section we analyze in details the general properties of the dynamical 
system such as transport and links to Lyapunov exponents. The second section is devoted to the construction of 
a thermodynamic formalism allowing us to define formal Gibbs measure, "partition function" and "free-energy" 
characterizing the stationary state. We discuss the thermodynamic limit in this setting, in the third section. We 
show in particular that, when the size of the system diverges, a Lee-Yang phenomenon [14] occurs, related to a 
loss of analyticity of the formal free energy. This unexpected result suggests that SOC systems might be closer as 
previously believed to classical critical phenomena and opens some effective way to "map self-organized criticality 
to criticality" [15]. We show that the Lee- Yang phenomenon is observed only when the energy is locally conserved 
(h = 0). Moreover, the critical exponent characterizing the divergence of the correlation length of the avalanche size 
distribution when /i ^ is analytically related to the angle that the Lee- Yang zeros do with the real axis with some 
analogy with usual critical phenomena. In section III B we also discuss the genericity of an important effect occuring 
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in the Zhang model: a weak form of initial conditions sensitivity. We wrote an appendix where a short presentation 
of the thermodynamic formalism, devoted to non-specialists, is presented. 

I. DYNAMICAL PROPERTIES OF THE FINITE-SIZE MODEL. 

A. Definitions. 

The Zhang's model is defined as follows. Let A be a d-dimensional box in 2**, taken as a square of edge length L 
for simplicity. Call N = #A = L'^, where # denotes the cardinality of a set. Each site i G A is characterized by its 
"energy" Xj, which is a non-negative real and finite number. Denote by X = {Xjj.g^ a configuration of energies. Let 
Ec be a real, positive number, called the critical energy, and M. = [0, Ec[^. A configuration X is stable when X e 
and unstable otherwise. In an unstable configuration the sites i such that Xi > Ec are called active. The dynamics 
action on X depends whether X is stable or unstable. 

If X is stable, one chooses a site i e A at random with probability j^, and add to it the energy ^ = 1 (excitation). 
If X is unstable, each active site loses a part of its energy, redistributed in equal parts to its 2d neighbors in the 
following way (relaxation). Fix two ^ real parameters, e G [0, 1[ and h > 0. Set 7 = e'*, e' = 67 = e^+'', a = ^^^d^ . 
When i is active it gives the energy aXi to its 2d neighbors and keeps the energy e'Xj. Therefore the energy is locally 
conserved in the case /i = whereas there is local dissipation ^ when /i > 0. If several nodes are simultaneously 
active, the local distribution rules are additively superposed, i.e. the time evolution of the system is synchronous. It 
is useful to write down the relaxation dynamics with the map: 

F(X) = X + (7 - l)eZ(X) * X + qA [Z(X) * X] (1) 

In this equation Z(X) is a A/' dimensional vector, such that Zi{X.) = if is stable and .^i(X) = 1 if is 
unstable. The * denotes the product component by component, that is if X, Y are N dimensional vectors, X * Y is 
the N dimensional vector of components XiYi. A is the discrete Laplacian on A. 

The succession of relaxations leading an unstable configuration to a stable one is called an avalanche. Let ^A be 
the boundary of A, namely the set of points in 'K'^ \ A at distance 1 from A. The sites of ^A have always zero energy 
(dissipation at the boundaries). Since the total energy in a finite lattice is finite, all avalanches stop within a finite 
number of iterations for /i > 0. The addition of energy is adiabatic. When an avalanche occurs, one waits until it 
stops before adding a new energy quantum. Further excitations eventually generate a new avalanche, but, because of 
the adiabatic rule, each new avalanche starts from only one active site. 

The structure of an avalanche is encoded by the sequence of active sites C = {Ct}i<t<T-c ^^^^^ = 
{j e A.\Xj > Ec in the t th step of avalanche} and where re, called the avalanche duration, is the smallest number 
such that Crc+i = 0- The step t = 1 corresponds to the initial excitation of some site i. Consequently either Ci = 



'^Note that the original Zhang's model corresponds to the case e = 0. The straightforward extension proposed here allows us 
to avoid pathological dynamical effects due to the existence of zero eigenvalues for the tangent maps when e = [7] . 

^The case h < would correspond to local energy injection. However, in this case there may not exist a stationary regime 
(see section IC). 
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(the addition of an energy quantum 6 is not sufficient to make the site i active) or Ci = {i}. Therefore, each avalanche 
can be labeled by a double index where the first index refers to the site where the energy is dropped and the 

second one to an enumeration of the possible avalanches starting at i (including the "empty" avalanche where the 
excitation of i does not render it active). The total energy of a stable configuration in a finite lattice being finite, 
(it is bounded by L'^Ec), the total number of possible avalanches is finite for finite L, (but diverges as i — > oo when 
/i = as discussed below). Call < oo the number of possible avalanches starting at i (note that depends 
on L and h). The size s of an avalanche is the total number of active sites, s{C) = ij^l^t Ct- The area a is the 
number of distinct active sites in C. Furthermore, one can use a, s, t to make a partial ordering of the avalanches : 
C ~<C' ^ n{C) < n{C'). The size s, area a, and duration r are examples of avalanches observables. We call avalanches 
observables a function n which assign to an avalanche a positive real number and such that C C C => n{C) < n{C'). 

Because all avalanches are finite (for finite L), and since we are not interested in the transients, one can, without 
loss of generality take all initial energy configurations X e A^. All trajectories starting from M. belong to a compact 
set v. Call M. = 'D\M.. M. contains the set of all unstable energy configurations achievable in an avalanche starting 
from an energy configuration m Ai. It is useful to encode the dynamics of excitation in the following way. Let be 

the set of right infinite sequences a = {a\, . . . , afc, . . .} , afe € A, and a be the left shif^ over S^, namely aa = 0203 

The elements of are called excitation sequences. The set x "D is the phase space of the Zhang's model and 
X = (a, X) is a point in the phase space. The Zhang's model dynamics is given by a map F : x X> — > x T> such 
that: 

XGA^=^F(X) = (aa,X + eaJ (2) 
XeM^ F(X) = (a, F(X)) (3) 

where eo is the a-th canonical basis vector of R^. The knowledge of an initial energy configuration X and of an 
excitation sequence a fully determines the evolution. One can endow Ej with a probability distribution correspond- 
ing to a random choice of the excited sites. The excited sites are chosen at random and independently with uniform 
probability. This corresponds to endowing E^ with the uniform Bernoulli measure, denoted by ul in the sequel. 

The model definition entails the convergence of the dynamics to a stationary state where the average incoming 
energy flux is compensated by the dissipated energy flux. It is furthermore conjectured in the SOC literature that this 
state is reached independently of the initial state. In the dynamical systems terminology this amounts to assuming 
that there exists a unique natural or Sinai-Ruelle-Bowen (SRB) measure [10-12]. Indeed, in this case, one has the 
following properties. Let fiL be an invariant measure on O. Let ip be some function, (integrable with respect to fiL)- 
Denote hy '4>l the time average : 



^Throughout this paper we will often think of the left BernoulH shift on as represented by the system z — > Nz mod 1, « € 
[0, 1]. This allows in particular to define the tangent map in the direction of the shift. 
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V;l(X)1^' lim i^^(F*(X)) (4) 

and let / il){'X.)dilL{X.) be the ensemble average with respect to fiL- Then /tj, is a SRB measure if 'ipLpQ = 
J ip{'X.)dfiLO^) fo'^ ^ of initial conditions X of positive Lebesgue measure [16]. If the SRB measure is unique 
then the time average and the ensemble average with respect to fiL are equal for a set of initial conditions of full 
Lebesgue measure. 

Equivalently, in this case, the SRB measure is the weak limit : 

^^L=lim^^j2^*\^yLX^,^°^) (5) 

°° t=i 

for any absolutely continuous measure fxf^^ characterizing the distribution of the initial energy configurations, where 
'F**{vl X M^^) is the image of vl x Ml ^ under the map F*. Note furthermore that if F is topologically mixing then : 

fiL = lim f\uL X nf ) (6) 

t — >CX) 

The existence of a unique SRB measure is somehow the minimal property to assume since it implies that, if one 
selects the initial conditions at random with a uniform probability (or any probability having a density), then the 
time average will not depend on the initial condition. Furthermore, the property (6) corresponds exactly to what is 
numerically observed. Starting from an arbitrary (absolutely continuous) probability measure /U^^ on the set of initial 
energy configurations and iterating the dynamics, the state F*(!/£ x /x^'), characterizing the statistical distributions 
of points at time t, converges to the stationary state jiL- This leads to state what we call the first SOC conjecture in 
the sequel: 

Conjecture 1 For any finite L, their exists a unique SRB measure [il obeying (6) for generic values of the parameters 
Ec,e,h>0 and whatever the lattice dimension d < oo. 

Note that the projection of /ii, on is the Bernoulli measure ul by construction. 

Proving this property in the Zhang's model is clearly a difficult task which is beyond the scope of this paper. We 
note however that this point has been widely discussed in a previous paper [7], where strong mathematical arguments 
in favor of this were given. Actually, the existence of a unique SRB measure was proved, but restricted to the one 
dimensional model, and to some Ec interval. Note a contrario that the failure of this property would lead to the 
coexistence of several stationary states depending on initial conditions. This has not been observed. 

Since the SRB measure is the stationary state, the physically relevant quantities characterizing the stationary 
regime are obtained from /ii. Call u)l '= Ai(^A ^ ■^)- By definition this is the probability to inject energy in the 
system. The average incoming energy flux, at stationarity, is given by the vector $i with components: 

= ^ = <t>L (7) 



*For any open sets Oi, O2 C A, 3T = T{Oi, O2) > such that Oi n F*02 #0, Vt > T [17] 
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corresponding to the average energy received by the site i, per unit time. Note that in this paper we consider, for 
simplicity, the case where the energy is uniformly distributed in the lattice, though most of the results hold in the 
more general case where is not spatially uniform. The quantity 4>l = X^iLi ^L{i) is called the injection rate. 
This is the average energy injected per site and per unit time. 
Define 

Pl{i) = Prob[Xi{t) > E^] = fiL [{X, Xi > E^}] (8) 

the probability that the site i is critical in the stationary regime, and let pL be the N dimensional vector, pL = 
The quantity 



av 

PL N 

i=l 



is often called the density of active site in the literature. It has been introduced by Vespignani et al. [18] as an order 
parameter in SOC models. Moreover, these authors introduced the quantity 

,,sdeldpl'' 

XL{h) = _ (10) 

that they interpret as a response function with respect to variations in the injection rate. In particular, considering 
the excitation at a given point and at a given time as a perturbation in the relaxation dynamics, xl (h) characterizes 
the linear response to this perturbation. We show below that, in our model, XL{h) diverges when h = and i — > oo, 
as expected for a critical phenomenon. 



In the conservative case h = 0, the numerical simulations report the following behavior [9]. Fix some avalanche 
observable n. Call PL{n) the probability distribution of n in the stationary state, for a box A of size L*^. It is observed 
that Pl (n) has a power law shape over a finite range, with a cut-off corresponding to finite size effects. As L increases 
the power law range increases. In the SOC literature it is assumed that the limit L — > oo of the model is well defined, 
and that the corresponding probability distribution for n, P*{n), is a power law. This is the ^ second SOC conjecture: 

Conjecture 2 As L ^ oo, for h = 0, for each observable n, Phin) converges to a power law 

P in) = , n = nn . . . 00 

^ ^ no- 
where no is model dependent. This entails the scale invariance of the corresponding stationary state which has 
therefore some common features with a critical state. r„ is called the critical exponent of the observable n. It is 
commonly admitted in the SOC community that a classification of the models can be made by the knowledge of their 
critical exponents ("universality classes"). 



^Note that the first (existence of a unique stationary state) and second (convergence to a limit state with power law statistics, 
as L — » oo) SOC conjecture are not inherent to the Zhang model but are implicit in all SOC models. The mathematical proof 
of the existence of a unique stationary state for a finite lattice has only been done in a few cases (such as the Abelian sandpile) 
[19], while the proof of the second SOC conjecture remains to be done for all the main SOC models. 
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When h^O, the energy balance changes since one has an additional dissipation (resp. injection) term when h> 
(resp. h <0). We show in the section ID that, when h> 0, the stationary state cannot converge to a critical state in 
the limit i — > oo. This regime is called subcritical. When h <0 one easily shows that there exists an h value, Hl <0, 
such that a stationary regime exists only if h > Hl- The corresponding regime is called supercritical. It exists only 
for finite L since /ii,— >OasL— >oo (see section I C). 

B. Transport properties. 

Wc now discuss the properties of the Lyapunov exponents and their connection with the energy injection and 
transport in the Zhang's model. This section presents original results which are an extension of the theory developed 
in [8]. Recall that the Lyapunov exponents are the logarithm of the eigenvalues of the + 1 x iV + 1 matrix 



lim. 



where " denotes the transpose. Note that the Osclcdcc multiplicative crgodic theorem [20] 
insures that the limit exists almost-surcly and is a constant. The particular structure of the dynamics (2) implies 
that DF-^ is block diagonal with a N dimensional block, DFx. corresponding to the relaxation dynamics and a one 
dimensional block corresponding to the excitation dynamics (sec footnote 3). 
The excitation dynamics is expansive and has a positive Lyapunov exponent : 

Xl{0) = ujLlog{N) (11) 

Note that Ai(0) is also the Kolmogorov-Sinai entropy. 

The Jacobian matrix D-Fx of F characterizes the energy transport during the avalanches, since the entry DF^ 
is the ratio of energy flowing from site j to site i, in t times steps, for the initial condition X [8]. Moreover, 
the relaxation dynamics is contracting [7] and, consequently, the N corresponding Lyapunov exponents arc strictly 
negative. Order them such that > Ai(f) > A/, (2) > ... > Xl{N). This hierarchy of Lyapunov exponents 
corresponds to a hierarchy of time scales for the energy transport in the lattice. Define M(X, i) = DF^.DF*^ 
and A — \imt^ oo M('X.,t)2t . M(X.,t) being symmetric it admits an orthogonal basis {vi(X,t)}^-^ and eigenvalues 
/ii(X, i). One can show that Xiii) = linii^oc 5(^09(/ij(X, f)). Furthermore, each Vi(X,t) converges exponentially to 
a vector Vi(X) in IR^, depending on X [21]. We call these vectors the Oseledec modes for the trajectory of X. They 
can be numerically obtained from the QR decomposition used in the computation of the Lyapunov spectrum (see 
[22]). The collection of Oseledec modes defines a hierarchy of nested subspaces (Oseledec splitting): 

= Vi(X) D V2(X) D...D Vjv(X) 

depending on X, where Vfe(X) = span |vfe(X), . . . , Viv(X)|, and such that the norm of a perturbation v G Vi(X) \ 
Vi+i(X) is given by : 

||£>F*5;.v|| =C(X,t)e^^W*||v|| (12) 

where limt^oo ^logC{X., t) = almost surely. Ai(z) is therefore the exponential rate of variation of ||v||. This expo- 
nent is closely related to the energy dissipation rate (see below and in [8]). 
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It has been shown in [8] that the spectrum of negative Lyapunov exponents can be roughly spHt into two parts: 
one part (fast modes and largest Lyapunov exponents in absolute value) corresponds to stabilizing modes contracting 
fast the dynamics onto the attractor; the other part (slow modes) corresponds to transport modes (see Fig. 1). The 
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slowest transport modes can be computed from a linear operator £ : R 

£(v) = V + 2(7 — l)epL * V + 2aA {pL * v) 



IR^ which acts on a vector v e as: 



(13) 



The Lyapunov exponents are given by the singular values of C. We have drawn fig. 1 the Lyapunov spectrum 
for L = 20, Ec = 2.2,6 = 0.1, h = 0.1, d = 2 and the singular values of (13), with (Fig. la) and without (Fig. lb) 
boundaries (in this case A is a 2 dimensional torus). One checks that the slowest modes are well approximated. 



Lyapunov spectrum ; L=20,Ec=2.2,E=0.1,h=0.1. 



Lyapunov spectrum (torus) : L=20,Ec=2,2,£=0.1,h=0.1. 
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FIG. 1. a. Lyapunov spectrum for L = 20, Ec = 2.2, e = 0.1, h = 0.1, d = 2 and the singular values of (13). Figure 1 b. The 
Lyapunov spectrum for the same parameters values, but without boundaries (torus). 



C can be interpreted as the transfer operator of a diffusion in a metric given by p^, with a probability 2(7 — l)epL 
of trapping a particle ®. In the presence of boundaries pL is non uniform in the lattice (see fig. 3a). The consequence 
is that the Lyapunov spectrum departs from a diffusion in a flat metric, except for the first mode (see Fig. la). A 
contrario, the same model defined on a torus with bulk dissipation has a spatially uniform p^. Consequently, the 
Lyapunov spectrum is similar to the normal diffusion in a flat metric, for characteristic times scales a bit larger 
than the average avalanche size. The Lyapunov exponents corresponding to smaller times scales departs from normal 
diffusion, due to the presence of a threshold in the dynamics definition (see fig. 1 b). 

This analogy with anomalous diffusive transport is fruitful since it allows us to interpret the Lyapunov exponents and 
the corresponding transport modes as the analogous of Fourier modes. In particular, the largest negative Lyapunov, 
Al(1), defines a characteristic time scale iL(l) — oi' microscopic escape rate, corresponding to the time taken by 

a particle to exit the system, either by disappearing in the bulk, or by escaping from the boundaries. In terms of the 
dynamical system (2) the largest negative Lyapunov exponent is the exponential decay rate of a generic perturbation 
about some point X. In particular the excitation of the site i corresponds to a perturbation about X oriented in the 



The factor 2 is due to the fact that a site cannot relax two successive time steps (see [7,8] for details). 
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i th canonical direction in IR" and belongs generically to Vi(X) \ V2(X). Consequently, Al(1) gives the decay rate 
of the initial local perturbation induced by the excitation and ^^(1) defines the characteristic time required for this 
perturbation to vanish. 

The first Lyapunov exponent is well approximated by (Fig. 2) : 

A.(l)--^ (14) 
and scales therefore like the dissipated energy per unit time (see theorem 2 in [8]). 



□ h=0. 
»h=0.1 




FIG. 2. First Lyapunov exponent versus L for Ec — 2.2, e = 0.1, h — 0,h = 0.1. In full lines are drawn the theoretical 

values given by (14) 



C. Stationarity equations. 

It is possible to obtain an approximate analytic expression for p^. It can be obtained by noting that : 

Fl - Xi - J {(7 - l)eZ(X) * X + aA [Z(X) * X]} d/iL(X) (15) 

is the average energy lost by each site within one step of the dynamics, in the stationary regime. At stationarity, this 
loss is compensated by the average incoming flux The balance equation writes: 

(7 - l)eU + aAU = -l>L (16) 

where IJ '= pL * Xj. The vector Xj = {X^{i)^'\'_^ is such that the i th component, X^{i) E [Xi\Xi > Ec]^ is 
the conditional expectation of Xi given that Xi > Ec- (Alternatively, X^{i) is the average energy of Xi given that 
Xi is active). 

It is easy to solve (15) by decomposing U on the eigenmodes of the Laplacian. pL can then be obtained by noting 
that the spatial fluctuations of Xj are small (except near the boundary). One therefore considers Xj as spatially 
constant. This yields the following "flrst order" approximation for '■ 

d 

Pl(x) = ^ A„ JJsm(fc^a;i) (17) 
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where x = [xi, . . . , Xn) is the set of coordinates of a point in a d dimensional lattice, n = (ni, . . . , n^) is the set of 
quantum numbers parameterizing the eigenmodes of the discrete Laplace operator, and 

d 

{-f -l)e + 2aC^cos{ki) - d) (18) 

is the corresponding eigenvalue with ki — 
The coefficients are given by: 



where E^^t = E^i^^(*) and. 



where n,- = 2to,- + 1. 



2: = ! 



The density of active site p°j^ — J2i=i is given by 



Pl 



-IL 



where: 



7L 



2 Y^ IlUcl 

L+l) ^ Sr, 



(19) 



(20) 



(21) 



The equations (17) and (20) are in very good agreement with the empirical values. We plot Fig. 3a the empirical 
and theoretical for L = 40, Ec = 2.2, e = 0.1, /i = 0.1 and Fig. 3b the empirical and theoretical curve p'^ as a 
function of L, for Ec = 2.2, e = 0.1 /i = and /i = 0.1. 



empirical + 




FIG. 3. a. Empirical and theoretical pL for L — 40, Ec — 2.2, e = 0.1, h = 0.1. Fig. 3b the empirical and theoretical (full 
lines) curve pj," as a function of L, for Ec = 2.2, e = 0.1 h = and h = 0.1. 
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The equation (20) is a balance equation. It can indeed be rewritten: 

Pl^l = h (22) 

where : 

corresponds (from eq. 22) to the average energy dissipated per unit time and per active site. It is called the dissipation 
rate in the literature [18]. 

Equation (22) implies that the susceptibility (10) obeys: 

.deL 



1 (pL dCL 1 

XL{h) 



d(j)L 



(24) 



ei, el d(j)L e-L 

Were the excitation rate per site ^i, and the dissipation rate ei, independent, would this relation reduce to 
= Yi- This corresponds to the result found by Vespignani et al. in a sandpile model where the excita- 
tion rate per site and dissipation rate were independent tunable parameters [18]. The relation (24) is an extension of 
this result to the case where Slj^l are not independent and tunable parameters but are determined by the microscopic 
evolution. (This is exactly what is meant by "self-organized"). 

As i — > 00 it is easy to see that the leading contribution in (21) corresponds to modes n = (m, . . . ,nd) such 
that Uk < P{L + 1) for some (3 > that can be taken arbitrary small. Indeed, for these modes, according to eq. 
(19), (25) Cn, ~ while s„ ~ (7 - l)e — (L+i)^ Sf=i '^i • Consequently, the sum on these modes scales like 
{L+lf'^ Y,* \ — rr where Y.* denotes the sum on the n's such that sup^. rik < /3{L + 1). On the other 



hand, one can decompose the remaining terms in the sum defining 7l into sums on the modes n such that exactly r 
components are smaller than /3(L + 1). Each sum scales like {L + 1)^'' where r < d. Therefore, as L — > oo, |7l| scales 
like fl 2^x-7)e h > 0, and like c{L + 1)'^^'^ when h ~ 0, where a, c are some non negative constants. Furthermore, 

E^g^ ^ L'^ since Ec < X^{i) < K, i — 1 . . . N, where K is some constant independent of L. Therefore, when L — > oo 
the dissipation rate. El, obeys the scaling relation : 

eL^X+[2{l-j)eA + CL-^] (25) 

where A,C are some positive constants. Consequently, when h > 0, the dissipation rate converges to a positive value 
as L —> 00, while it converges to like when /i = In Fig. 4a we plotted the dissipation rate cl- In full 
lines are drawn the fitting curves X^CL~^ (conservative case h = 0) and [a + CL~^] (non conservative case 
h = 0.1) corresponding to eq. (25), where A,C have been obtained by fitting. We found a = 0.0238 ± 0.0001 and 
C = 5.959 ± 0.14. In Fig. 4b ujl is plotted . 
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FIG. 4. a Dissipation rate versus L for Ec — 2.2, e — 0.1 in the case h — 0, h — 0.1 (insert; corresponding values of X^). In 
full lines are drawn the fitting curves X^CL~^ (conservative case ft = 0) and [a + CL~^'^ (non conservative case h = 0.1) 
corresponding to eq. (25), where A,C have been obtained by fitting. Fig. 4b The corresponding graph of ljl- 



Note that there exists a stationary regime only if the dissipation rate is positive. Henceforth, there exists a Hl given 
y En sL = such that 
and converges to as i — > cxd. 



by En ' = such that, for h < there is no stationary regime, behaves asymptotically like — ^j^^^-^^ — - 



D. Avalanches observable distributions. 

Call PL{n, h) the probability distribution of an avalanche observable n, in the stationary regime, for a lattice of size 
L. Call 

<n>L,h^ nPL{n,h) (26) 

n=0 

the average value of n. ^£(^) is the maximal value that the random variable n can take, among all the avalanches 
having a non zero probability of occurrence in the stationary regime, in a lattice of size L. Since the energy used to 
initiate an avalanche is dropped locally, C£(/i) is the largest scale where the effect of a local perturbation is observed, 
within one avalanche. In this way it corresponds to a correlation length within one avalanche. This is a function of 
h. As discussed above < oo for > 0, L < oo. To simplify the notation we will set : P]^{n,0) = Phin) (resp. 

<n>L.o=<n>L, Cm=il). 

In the subcritical regime /i > the avalanche size (resp. duration, area, etc..) is bounded VL. Indeed, consider the 
evolution, under the singular diffusion (1), of a configuration Y in the infinite lattice^ Z'', such that Ec > Yi > Ec — i], 
?7 > 0, Vi G Z'', with excitation at some point ia G 2''. Consider simultaneously the normal damped diffusion 
Y'{t + 1) = Y'(<) + (7 - l)eY'{t) + aA [Y'{t)] on Z"* with a source 6^1 applied at time t = at the site io, where 
y/(i) ^ maxiY.it) - {Ec - rj), 0). By definition Y,{t) < Y^'{t) + Ec - 'n,Vt > 0. Y'{t) converges asymptotically to 



''With bulk dissipation since h > 
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and it is easy to show^ that there exists a bounded region TZ C 2'' containing iq, and a time to such that, Mt > to and 
Vi G Z''\7^, < l^'(t) < f . Since Yi{t) < Yl{t)+E^-'q, Yi{t) < E^-^ < E^,yt> to, Vi €71^X11. Consequently, the 
avalanche C(Y) does not got beyond TZ and the duration (resp. area, size) of C(Y) is bounded. Now, since the largest 
avalanche size in a finite box A c 2'' is generated from a stable configuration X such that Xi < Ec it is possible to 
find a configuration Y as above such that C(X) c C(Y), provided r] is sufiiciently small (take t] < Ec — maxi^\Xi). 
Consequently, < n(C(Y)) < oo,Vi,V/i > 0. It follows that the maximal size, area, etc... are bounded together 

with the corresponding average values when h > 0. Moreover, the energy excitation rate, u!L{h), scales like ^^^^ ^ - 
Consequently, u)L{h) converges to a strictly positive value in the subcritical regime (see Fig 4b). 

On the other hand, the average avalanche observable and the corresponding correlation lengths diverge in the 
conservative case, when L — > oo. Firstly, the behavior of the average avalanche size < s >L,h sis a. function of h and 
L is obtained from a stationarity condition insuring the balance between the incoming energy flux (5 = 1 is injected 
each time an avalanche ends) and the average energy flux dissipated within an avalanche (< s >L,h e.L{h)). This 
gives: 

< ' ^''"^ ilW ^ X+[2(1-7U + CL-] ^''^ 

Therefore, < s >L,h converges to a constant when h > 0, while it diverges like in the conservative case [8]. The 
adiabaticity condition imposes that the number of active site at the t-th avalanche step is bounded from above^ by 
X^^^i(2rf)*^~^. This implies that the expectation of the avalanche observables a, r also diverges. It is observed that 
<n >L diverges like: 

<n>L-^L''- (28) 

where n = a,s, t, and 7„ > in the conservative case. 

Since most of the results available (in particular the numerical ones) are obtained for finite L, one has to propose a 
finite size scaling form allowing to extrapolate to L — > oo from finite size lattices results. The most common scaling 
form (see in particular [18]) is the following (adapted to our notations): 

Tl 

PL{n, h) = n-^"/(^) , n > no (29) 

where no > is model dependent. 

This scaling form is discussed in more details in section II D, but this is obviously the analog of the finite size scaling 
forms used in critical phenomena. The form (29) implies < n >l,/(~ CEC^)^^"^"^- Therefore, in the conservative case, 
C£(0) ~ L'^", where: 

^" 2-r„ 



*On TL"^ the eigenmodes of the damped diffusion equation are X{k) = (7 — l)e — ak^ where k is the wave vector in the Fourier 
space. The presence of the damping coefficient (7 — l)e insures the existence of the region IZ. 
^This is certainly not an optimal bound since, in particular, a site cannot relax two successive time steps. 
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and Phin) = "^^i^ corresponds to the Kadanoff et al. [23] finite size scaling form (61). Clearly, in this 

scaling form, t„, (3n are characteristics exponents allowing to determine the universality class of P*{n). 

The form (29) also implies that the correlation length of avalanche size is a function of the dissipation rate with a 
scaling form : 

a{h) = eL{hr^ (30) 
where a = 2 — Tg. Therefore, for h = 0, the correlation length diverges like L^" where (3s = > ^^d, for h > 0, 

— — A 1 

^£(/i), converges to a constant ^^{h) = (2A(1 — 7)e) " . Finally, as /i ^ 0, ^^{h) ~ h~°, where S = ^. Consequently, 
^ is a characteristic exponent related to the singularity of correlation length ^£(/i) as /i — > 0. In the section III A we 
show how exponents r„, /3„, ct can be determined from the behavior of the zeros of a suitable partition function in full 
analogy with usual critical phenomena . 



Prom eq. (20,25) the density of active site behaves asymptotically like: 



Consequently, the density of active site converges to whatever h>0. It scales like — ^ ^ L ^ for /i > and like 
^ for /i = 0, where <t>l^L'^-. 

Finally, let us discuss the behavior of the susceptibility (24). It behaves like ^ =< s >l as i — > oo unless the 
term p'l^^ converges to 1 as -L ^ oo. This extends the relation already found by Vespignani at al. [18] for sandpiles, 
to the Zhang model. 



E. Extrapolation to L — » oo. 



In this section we want to make some remarks, based on the previous analysis, concerning the limit of the model when 
the size L tends to infinity ( "thermodynamic limit" ) . Note first that a mathematical definition of the thermodynamic 
limit raises serious problems in SOC. In statistical mechanics the finite volume equilibrium invariant^^ measure is 
known: this is the Gibbs measure ' where the finite volume Hamiltonian Hy y contains interactions with the 

exterior V = TL'^ \ V. In SOC models the stationary state is the result of a specific non Hamiltonian microscopic 
dynamics and its explicit form is not known. Moreover, in models like the Zhang's model, there are typically infinitely 
many ergodic measures and one has to add additional constraints in order to select the physically relevant one. 



"To avoid this case it is enough to assume that 



< oo, VL since 0. As noticed above cl, (pL are not independent 

parameters and are constrained by the dynamics. They depend on the stationary state which is itself determined by the control 
parameters Ec, e, 5. It is therefore sufficient to assume that the variations of eL,4'L are regular with respect to variations in the 
control parameters. This holds unless the system is at a bifurcation point. 

^^Note that the form of the equilibrium measure is in general obtained from general principles, without reference to the 
underlying microscopic dynamics, though dynamical systems admitting the Gibbs measure as the unique invariant measure 
can easily be constructed (Glauber dynamics for example). 
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In statistical mechanics the infinite volume Gibbs measure can be constructed via the Dobrushin-Landford-Ruelle 
(DLR) specification [24]. One attempts to find a measure g* on (or R'') such that the conditional measure Qyy, 
restricted to a finite volume V, with specific boundary conditions in V, is the corresponding finite volume Gibbs 
measure. Under suitable conditions on the interactions, one can show that this measure exists and is unique [25]. A 
similar construction for SOC models can be conceived but raises technical difficulties. One would like to construct 
a measure fi* on 2'' such that the conditional measure /i^ with zero boundary conditions on A = 2'' \ A is the 
stationary state of the finite volume model {fiL in our case). However, the construction of a SOC dynamical system on 
Z*^ (or a Markov process) and the computation of the corresponding invariant measure has only been done in specific 
one dimensional sandpile models [26] , and has not been yet done in larger dimension for the main SOC models (BTW, 
Zhang, OFC, etc.). The main obstacle in the construction is the presence of long range correlations which prevents 
the use of classical theorems [26] . 

In the present paper, we will therefore not try to construct the thermodynamic limit in a rigorous way. Rather, we 
will attempt to build a finite size thermodynamic formalism for the Zhang's model with the objectives presented in 
the introduction. Note that when dealing with the notion of thermodynamic limit one is usually faced to the delicate 
question of the order in which the limits L — > oo, T ^ oo are done. However, since usually the scaling properties 
of the stationary state w.r.t. L are discussed in the SOC literature, the limit i ^ oo is implicitly done after the 
limit T ^ 00. This is what we will do in this paper. Note that there are no a priori reason why these two limits 
should commute (usually, in critical phenomena, they do not). Note that the case T^oo,L^oo,h = 0, there is no 
dissipation at all in the limit model and consequently no stationary regime. 

Let us now make several important remarks coming out from the analysis made in the previous sections. For 
h ~ 0, when L ^ oo the system reaches a state where the correlation lengths = ^) ''') ^ diverge, and where the 

avalanche are statistically distributed according to a power law. In this sense, it corresponds to a critical state, which 
is reached spontaneously by the only effect of the dynamics, the adiabaticity condition in the energy injection, and 
the vanishing of the boundary dissipation rate. From the dynamics point of view, the positive Lyapunov exponent 
u)L ^og{N), which is also the Kolmogorov-Sinai entropy, vanishes since the excitation rate u)l tends to zero. Therefore 
the asymptotic state has zero entropy. In the same time, the dissipation rate vanishes and correspondingly, the 
first negative Lyapunov exponent tends to zero (see eq. (14)). Consequently, the Zhang's model loses its hyperbolic 
structure in the limit L ^ oo. This is clearly expected since the loss of hyperbolicity is a necessary condition to 
have an algebraic correlation decay. Indeed, hyperbolic dynamical systems have exponential correlation decay. A 
sufficient condition is given by the vanishing of the spectral gap in the spectrum of the Perron-Frobenius operator 
(see the discussion). For h > 0, Al(1) still converges to but the positive Lyapunov exponent (the entropy) diverges 
like log(A^). On the other hand, the discussion above shows that there is no critical state in the thermodynamic 
limit. Hence, h can be used as a control parameter allowing to tune the system to the critical regime. In the next 
sections we show that the corresponding phase transition can be handled, in the proper thermodynamic formalism, 
with the analysis tools of statistical mechanics (Lee- Yang zeros). The next section is devoted to the construction of 
the thermodynamic formalism for the Zhang's model. 
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II. THERMODYNAMIC FORMALISM. 



The evolution equation (2) characterizes a deterministic dynamical system. One can also describe the dynamical 
evolution of the Zhang model by a Markov process with an uncountable phase space. However, the hyperbolicity of 
the finite size model, allows us to reduce the dynamics to a Markov chain with an most countable phase space. Indeed 
there are energy configurations equivalent from the dynamical point of view, in the sense that their trajectories are 
asymptotically indistinguishable. From the sets of equivalent configurations one can build partitions of the phase 
space allowing to encode the dynamics symbolically by an at most countable Markov chain. Moreover, in some cases, 
the Markov chain is finite. 

Furthermore, by redefining the dynamics in terms of return maps, we construct the symbolic dynamics in such a way 
that each symbol corresponds to energy configurations undergoing the same avalanche when a specified site is excited. 
In this way, the coding is both relevant to characterize the microscopic evolution but also the avalanche dynamics. The 
first SOC conjecture implies that there is a unique invariant measure in the code space, for the corresponding Markov 
chain, corresponding to the SRB measure. By definition, it characterizes the microscopic dynamics. In particular its 
support corresponds to the so-called "SOC attractor" [1,3]. But, by construction, it also characterizes the macroscopic 
avalanches distribution. This measure is a particular example of Gibbs measures (see appendix). We construct other 
examples of Gibbs measiire in subsection II D, by changing the weights in the Markov transition matrix and discuss 
their connexions with the microscopic dynamics and the distribution of avalanche observables. 



We first redefine the Zhang's model dynamics in terms of return maps. The avalanche after the excitation of a site 
maps unstable to stable configurations. One can view this process as a mapping from M. ^ M. where one includes 
the process of excitation. Let Tj be the map A\ ^ M. which associates to a stable energy configuration X the next 
stable configuration resulting from an avalanche obtained by the excitation of the site i in the configuration X. Then: 



in terms of the the mapping (1), where r(i,X) = inf {< > 0,i^*(X + 5ej) € M.} is the duration of the avalanche 

obtained by exciting the site i in the stable configuration X. 

From the mappings Ti^i G A, we construct a new ("Poincare like") dynamical system where the phase space is 
= X Al, where, as above, is the set of right infinite excitation sequences a = {at, . . . , a^, . . . |afe e A} but 

X is now a stable energy configuration. X = (a, X) is now a point in Q.. Let tti and 7r2 be the canonical projection 

respectively on and M (namely 7ri(X) = a, 7r2(X) = X). 

The evolution is determined by a dynamical system of skew product type, T : ^ — > ^ such that: 



A. The return maps. 



T,(X) =i^^"(*'^)(X + <5e,), XeM 



(32) 



T(X)^=i'(aa,T„,(X)) ; 



X^=i'(a,X) 



(33) 



Set X(t) = T*(X) and X(f) = 7r2(T*(X)). X(t) is now the stable energy configuration obtained after t avalanches, 
starting from the energy configuration X when the excitation sequence is a = tti (X) . 
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Due to the particular structure of the mapping (1), each map Tj, i e A, is a piecewise affine map [7]. Denote by 
T(^i,j) the j th afRne component of the map Tj. Call Af(j the domain of definition of T(^i,j)- There is a one to one 
correspondence between the set of avalanches and the set of afhne mapping: T(^i,j) maps the energies configurations in 
M(i_j) to the stable configurations resulting from the avalanche (Recall that each avalanche can be labeled by a 

double index where the first index refers to the site where the avalanche starts and the second to an enumeration 
of the possible avalanches starting at i). For each i G A the M(ij)'s, j = 1 . . . constitute a partition of M. 

The set of points <S where T is not continuous is called the singularity set. This is an union of hyperplanes 
constituting the borders dM^ij^ of the Mj^ j)'s. In more mathematical terms: 

S=[j[jdM^,j^ (34) 

where: 

dM(^ij) = {X e A^, 3keA,3te{l... T{i, X)} , F^(X + ^e^) = E^} (35) 

The singularity set S is therefore the set of stable energies configurations such that some sites have an energy 
exactly equal to Ec at some time during some avalanche. These configurations have therefore some site "right at the 
threshold" at some time of their evohition. They are therefore highly sensitive to an infinitesimal perturbation on the 
energy of these sites. S plays an important role in the dynamics, discussed in the next section and in section IIIB. 



It can be proved that each map is a quasi contraction. More precisely it contracts on the subspace generated 
by the the canonical basis vectors corresponding to the active sites, and is neutral on the remaining part of IR^. 
Furthermore, 

det(T(,,,.)) = e««^'^-)) (36) 

where s{{i,j)) is the size of the corresponding avalanche [7]. 

Moreover, det{TTi{DT^)) = log(iV) since the shift on is conjugated to Nxmodl. Consequently, 7ri(T) is expansive 
with a positive Lyapunov exponent Ci(0) = ^og{N). It has also L'' negative Lyapunov exponents Cl(*)) i = 1 • • ■L'^. 
The Lyapunov exponents Cl(*) are related to the Lyapunov exponents of section IB by a simple reparametrization of 
time. Indeed, in the dynamical system (33) a unit of time corresponds to the time elapsed between two excitations. 
In the dynamical system (2) the average corresponding time is Jj- . Consequently, the Lyapunov exponents Cl (*) are 
given by: 

CL{i) = ^, i = Q..-N (37) 



^^The maps T^i,j) play in some sense the role of the toppling operators introduced by Dhar for Abelian sandpiles [19]. However, 
the structure of the Zhang model is more complex since these operator do not commute and since they do not preserve the 
Lebesgue measure. 

^^The number of singularity planes is finite when L < oo but it becomes infinite as L — » oo when h = 0. 
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Since CI has a product structure, and since the excitation is independent of the dynamics on energy configurations, 
the invariant measures we consider decomposes as fi^ x where fi^ is the induced measure on the unstable direction 
or excitation measure and /if, is the induced measure on ^4 or measure on the stable energy configurations. In the 
Zhang's model fi^ is the uniform Bernoulli measure since the successive excited sites are chosen independently with a 
fixed rate ^. To the SRB measure defined in section IB corresponds a SRB measure for (33), with product structure 
/t^ X /t£. In the sequel we will keep the notation fiL for the SRB measure of (33). 

B. Symbolic dynamics. 

In this section wc build a symbolic coding of the dynamics, relevant to characterizes both the microscopic dynamics 
and the avalanches dynamics. For, one has first to build a suitable partition V of the phase space = x A4. 

A rooted avalanche is a pair (C, B) where C is an avalanche and B d Q such that for all points X G _B the only 
possible avalanche is C. Equivalently, if C = {i,j) then tt2{B) C Af(i.j) and tti{B) = [i], where [i] denotes the 
set of excitation sequences whose first digit is i (cylinder set) |a e \ ai — i^. A system {C^, B^)^^j of rooted 
avalanches, where the set of indexes 2 is possibly countable, is called complete if [J^f^j B^^ — Cl and B,^ H B'^ = 
for Lo ^ Lo' . Hence, {Bu>}^^-x is a partition of CI. By construction to each element uj of this partition corresponds a 
unique avalanche [i,])- Equivalently, to each lj one can associate a double index (i {uj) ,i {lj)) where the first index 
refers to the site where the avalanche starts and the second to the corresponding avalanche. Since tt2{B) C M^ij-^ the 
avalanche corresponds in general to several w's. However, the partition Vnat = {[*] x -^(^j) I^ga j-i called 
the natural partition in the sequel, is such that there is a one to one correspondence between the symbols uj and the 
avalanches. By definition any complete system of rooted avalanche corresponds to a partition generated from Vnat- 

To each complete system of rooted avalanches (Cc^, Bi^)^^-j-, we associate a transition matrix A = {Aij)i_jex such 
that Aij = 1 if Bj n T~^{Bi) ^ (the transition i ^ j is, legal) and Aij = otherwise, and a transition graph 
Q with vertices lo £ I and oriented edges i —>■ j for all pair (i,j) such that Aij = 1. This provides a symbolic 
dynamics description of the Zhang's model where the trajectory of a point is represented by a legal sequence of 
symbols . . . luiuj2 ■ ■ ■ w„ . . . corresponding to the partition elements that this point meets in its history. However, 
except for special partitions called generating partitions, the coding is in general not unique. 

The transition graph Gnat of the natural partition characterizes the legal compositions between avalanches (resp. 
mappings). However, though the natural partition is the most suitable to describe the avalanches evolution it does 
not have, in general, the Markov partition property allowing to encode the dynamics in a Markov chain. We say that a 
complete system of rooted avalanches (resp. the corresponding partition {Buij^^^x) Markov partition property 

if, for all couple (i, j), such that Aij = 1 : 

(i) 7ri{T{Bi)) D ni{Bj) (38) 

(ii) TT2{T{Bi))dTZ2{Bj) (39) 

By construction the property (i) always holds. If (ii) holds then the corresponding complete system of rooted 
avalanches has the following properties: 
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• (i) Every (legal) finite path a;i,...,a;„, iVk S I corresponds to a legal sequence of avalanches 

j(a;i)), . . . , j(a;„)). 

• (ii) For any X each orbit segment < T' (X) [ is realized as a path of length n on the Markov transition 

I J l<l<n 

graph. 

Except for a non generic set, the trajectory of any point X in the phase space O is uniquely encoded by an infinite 
sequence of symbols corresponding to the partition elements visited by X. 

The physical interpretation of the property (ii) is that, starting from any energy configuration in one element of the 
partition and exciting a given site the next avalanche is uniquely determined by the next excited site. Consequently, 
the natural partition Vnat docs not obey (ii)) in general, though, in some specific examples, there exists a finite iterate 
of T, T* such that T^iVnat) has the Markov property. This is the case of the one dimensional Zhang model where 
Ec e]l,2],/i = 0,e = [7]. In the general case, in order to have a coding relevant both for the avalanche dynamics 
and having the Markov property, the strategy is to build a partition from the natural one, by cutting it in sufficiently 
small sub pieces. However the resulting partition can be countable. 

The condition (ii) is rather strong since it deals with all points in O. Since we essentially want to have a symbolic 
description of the asymptotic dynamics, namely on the attractor which is a subset in O, it is in fact enough to prove 
the existence of a finite Markov partition in the usual sense [17]. This is insured by the existence of a local stable 
manifold of sufficiently large diameter, for almost every point in [7]. Recall that the local stable manifold WfQg(X) 
of X G is the largest connected set such that VY G Wf^^, d(T*(Y), T*(X)) — > when t ^ oo, d being some distance 
on ri. Therefore, the trajectory of any point belonging to Wf^^pC) is asymptotically identical to the trajectory of X 
and these points are equivalent from the dynamical point of view. Equivalently, if X has a local stable manifold of 
diameter larger than some i] > then a small perturbation of size < t] will be asymptotically damped (see section 
IIIB for a physical interpretation of this effect.). 

Were the map T be regular, were the existence of local stable manifolds and Markov partition be insured by 
the standard theorems on regular (C^+") hyperbolic dynamical system [21]. However, T is not continuous on the 
singularity set S and some points have no local stable manifolds. The main problem is to estimate the fiL measure of 
the "bad" points having no local stable manifold. The following result is useful. 

Proposition 1 Let Us{S) — |y g U \ d{Y,S) < be the S -neighborhood of S. Assume that for 6 > sufficiently 
small: 

fiL [Us{S)] < C5" (40) 

some C > 0, some a > 0. Then (il almost every point has a local stable manifold of positive diameter. 

Proof The existence of a local stable manifold of positive diameter for a point X in the support of /ti, is ensured 
if one can find a number > 7 > Ai,(l), and a time to < 00 such that \/t > to, rf(T*X, <S) > e'*'*. Indeed, in this 
case a sufficiently small ball around X is asymptotically contracted faster than it approach the singularity set, and 
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consequently, all the points of this ball belong to the stable manifold of X. Consequently, the set of bad points is 
included in the set : 

C30 CXD 

{X I Vio > 0, 3i > to, d{T*X,S) < e-^*} = f| |J T"* {lie-,t{S)) 

to=0 t=to 

for some 7 such that > 7 > Ai:,(l). Prom the Borel-Cantelli lemma the measure of this set is zero provided the series 
YltLoi^L [Uen{S)] converges. This is true if the condition (40) holds. □ 

We have not established yet a mathematical proof of the condition (40) but we have strong numerical evidences 
presented in the section IIIB. In the sequel we will therefore assume that fiL almost every point as a stable manifold 
of positive diameter. Then, it is possible to generate a finite Markov partition V encoding the asymptotic dynamics 
and such that each symbol corresponds to a unique avalanche. Note however that several symbols may correspond 
to the same avalanche. Nevertheless, though an avalanche can be represented by several partition-elements, 
in order to simplify the notation we will however identify w and the avalanche {i{u)),j{w)) whenever it makes no 
confusion. Denote by a) = a;_„ . . . a;_ia;oa'ia;2 . . . a;„ . . ., where A^„,a;„+i = 1, n € Z, a legal bi-infinite sequence. Call 
X (resp. X~^) the set of bi-infinite (resp. right infinite) legal sequences. There is a one-to-one correspondence (up to 
a set of zero Lebesgue measure), denoted by ~, between a symbolic sequence w and a point X in the phase space. 
Let TT = (7r+,7r~) be the corresponding conjugating mapping, such that if a) ~ X = (a, X) then 7r+((D) = a and 
7r~((D) = X. Note that 7r+ projects on the unstable direction (S^^), while tt" projects on the stable space {M). Write 
[u]n for a n-cylinder (this is the subset of X where the sequences have the same n first digits as u). Denote by the 
shift on X. The forward orbit of u) under encodes the excitation sequence by definition and the backward orbit of 
u) encodes the point in the energy configuration space A4. 

C. Markov chain and SOC state. 

We construct a Markov chain by defining a transition kernel W from the matrix A. A has the following property. 
Let £)+ be the set of follower elements of u (and the set of predecessors of to). #1?^" > N, Vw since after the 
avalanche corresponding to lu one can excite any site in A. Moreover, the Markov partition property (38 ii) implies 
that if /3, 7 G and /3 7^ 7, then i{(3) ^ 2(7). Hence, #-Dtj — N, Vw. Consequently, A has exactly N non zero 
component per row, corresponding to each of the TV possible choices of an excited site. 

Since the excited site of A is chosen with respect to the uniform Bernoulli measure the transition probability of the 
Markov- Partition for any edge from u) to is constant, equal to j^. Therefore the transition kernel is : 

W=^A (41) 

As a consequence, W is a sparse matrix. Indeed, the number of symbols iv is at most equal to the number of 
possible avalanches, which is bounded from below by the largest avalanche size (resp. duration, area) Since, for 
/i = 0, ^£ ~ L^", ps > d, the proportion of non zero entry on each row, ri, = ^L'^~^'', tends to zero as L ^ 00. 
This remark justifies somehow to approximate the transition matrix by a random matrix. This could be used as an 
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approximation to determine the spectral gap of W (see the discussion). 

The set of symbols T decomposes into transient and recurrent nodes. Call TZ{T) the set of recurrent nodes and 
-A/z, '== #TZ{I) the number of elements in 'R{I). The first SOC conjecture requires that the set of recurrent node is 
irreducible. This is not in contradiction with the sparseness of W. For example, it can be proved that for sparse 
random matrices, with K{m) uniformly distributed non zero entry per row where m is the size of the matrix, the set of 
nodes is almost-surely constituted by one irreducible recurrent cluster as m tends to infinity provided K{m) > log(m) 
[27]. 

The equivalent of the convergence property (6) is : 

mL= lim vW* (42) 

where v is any initial probability distribution on T. This property holds when W is mixing In this case W has a 
unique eigenvalue 1 and a unique left eigenvector tul, corresponding to the invariant probability distribution of the 
Markov chain [28]. mi,(a;) is non zero only if w e TZ{1). Furthermore, there is a spectral gap between the second 
eigenvalue with the largest modulus and 1. The gap gives the exponential correlation decay and also the rate of 
convergence to equilibrium. It is expected that the gap vanishes as L — > oo for /i = and stays positive for h> 0. 

The invariant probability distribution tul characterizes the probability of occurrence of any recurrent symbol at 
stationarity. This is therefore the fundamental object characterizing the SOC state. In particular the probability 
distribution of an avalanche observable is given by : 

PL(n,/i)= ^ niLico) (43) 

uj^X ,n{uj)=n 

Recall that n{w) stands for n{i{u)),j{w)). This is the value that the observable n takes in the avalanche i{u)),j{u)). 
Note that mLjfiL depend on h but we dropped this dependence in order to simplify the notations. The moments of 
n are given by: 

tJGl n=0 wel,n(w)=n "=0 

More generally, the joint probability /xl on the space of infinite symbolic sequences X is obtained from the Chapman- 
Kolmogorov equation. For any cylinder [u)]t = i^i ■ ■ ■ ojt- 

Hiii^W) = mL(wi)>Va;iW2>Va,2a>3 • • • W^^t-H^t (45) 

Note that fiL is the measure of maximal entropy \og{N). 

The measure /ii projects down to the SRB measure fiL on via fiL = fiL°T^~^ with marginals = /ii o 7r+ and 
/*£ = p,LOTT~ . For the following we will therefore make no distinction between the average with respect to jjL or with 
respect to fiL. From the measure fiL one can compute all the time correlations (therefore the joint probability on the 



There exists n > 1 s.t. for each pair i,j W"{i,j) > 0, i.e. there is a path connecting each node 
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space of trajectories of n{t)) of the observable n (see eq. (56)). 

The knowledge of the transition matrix W allows us to determine the invariant measure ttil and the spectral gap. 
However its structure is not known in general. An interesting, but very specific case, occurs when the number of 
predecessors of w, #-Dj , is also equal to N, Vw. This happens for example in the one dimensional Zhang model where 
Ec 2], e = 0, /i = [7]. Then W is a bi stochastic matrix and mi, is the uniform measure, mL{co) = j^, on TZ{T). 
However, in general, #D~ depends on w. 

D. Thermodynamic formalism, Gibbs measures and generating function of avalanche size distribution. 

In this section wc use the thermodynamic formalism (see appendix) to construct Gibbs measure, by choosing 
different families of potentials. These measures arc relevant to characterize microscopic properties (like the fractal 
structure of the SOC attractor), but also the avalanche distributions. 

The SRB measure /t^ is the equilibrium state which maximizes the potential - log(det(7r+(i:)7^))) [16]. In the 
Zhang model this potential is a constant — log(A^). It follows that the SRB measure maximizes the entropy. Since the 
maximal metric entropy is the topological entropy log(iV), it also follows that the SRB measure has zero pressure. 

Due to the product structure of the Zhang's model: 

log {det{D%)) = log(7V) + log {det{DT^,)) (46) 

where DT^^^ stands for DT^^^^^^^j^^^^-^y Let us introduce a first family of potentials: 

(l>l{Lb) = - \og{N) + (3\og {det{DT^,)) (47) 

and denote the corresponding family of equilibrium states by up in the sequel. The potential (j)^ are Bernoulli 
potentials (they depends only on the first symbol in cD). Consequently, they trivially match the condition (72) in the 
appendix. 

The case /3 = corresponds to the SRB measure. The quantity \og{det{DTi^-^)) is the volume contraction in the 
space of energy configuration, induced by the relaxation dynamics. It is related to the energy transport in the lattice 
and to the avalanche size s{uj\). Indeed, eq.(36) implies: 

4(d;) = - log(7V) + /?log(e)s(wi) (48) 

The corresponding partition function (73) writes: 

Zi(/3)=7V-^ ^ e''Sr=i^(-*) (49) 

The corresponding topological pressure is: 

rm"^ Y,rn ;^log(Zi(/3)) (50) 

Then from eq. (76), 
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d(3 



(51) 



3=0 



This equation is analogous to the one giving the magnetization when deriving the free energy w.r.t. to a uniform 
local field. It relates the average volume contraction rate given by the sum of negative Lyapunov exponents Ci(*)> * = 
1 . . .L'^, to the average avalanche size [7] 



N 



(52) 



i=l 



The second derivative of the free energy (50) with respect to /? is: 



(3=0 

This oqiiation gives therefore the variance, along a typical trajectory, of the Gaussian fluctuations of the avalanche 
size (rcsp. contraction rate) around the average values. In section III A we discuss the failure of differentiability of 
for /i = 0, in the thermodynamic limit, which implies that (53) diverges and that the central limit theorem is 
violated. 



More detailed informations about the probability distribution of avalanches size, or more generally of any other 
avalanche observable n, can be obtained by introducing in the potential the formal equivalent of a local magnetic field 
rj which allows, in statistical mechanics, to compute the local magnetization and higher order cumulants. A natural 
extension of (47) is therefore the time dependent potential: 



(j)n{t, cD) = - \og{N) + r]tn{uit) 
where rj = {rji, . . . ,rit, ■ ■ ■), and the partition function : 

This expression reduces to (49) when rjt = /31og(e), n = s. Set : 

^Liri) = lim ^ log(ZT(?7i, . . . , tjt)) 
Thirj) generates the fc-time correlations of the observable n: 



drji... drjk 



(ni;n2; . . • . 



t;=0 



In particular, the cumulants CL{k) of the marginal distribution Phin) are given by: 



dr]f 



= CL{k) = 



r)i=0 



dt'' 



t=0 



where : 
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GL{t) log((e*")^ _,) = log( J2 PL{n, h)e'n 



(54) 



(55) 



(56) 



(57) 



(58) 



(59) 



n=0 
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is the generating function of the cumulants of Phin, h). 

Thirf) plays a similar role as the generating functional in quantum field theory and r] acts as a conjugated field 
to the observable n. It also corresponds to a time dependent perturbation of the trajectories or, equivalently, to 
a time dependent potential. It must in particular be emphasized that the topological pressure in (56) contains all 
informations about the scaling property of the probability distribution on the space of trajectories of the avalanche 
observable n. Consequently, a scaling theory, analogous to statistical mechanics seems achievable. Formally, when 
h = 0, one can associate to each local field r]i an exponent yi and seek a scaling form for J^Liv) such as : 

J^l{v) = LP-Hiiy^rn, . . . . . .) (60) 

where W is a universal function. This provides the scaling of the n points cumulants. However, one has introduced 
an infinite number of fields and this formula raises the question: what are the relevant fields (in the renormalization 
group sense, that is unstable directions for the renormalization flow) or in other words which scaling exponents 
define the universality classes in SOC ? In classical critical phenomena the renormalization group approach, and the 
renormalization properties of the Hamiltonian and free energy [29] allow to select the relevant fields and two scaling 
exponents are extracted; a, the specific heat exponent and rj (correlation length). In the SOC case the main difficulty 
is to adapt the renormalization schemes and to select the relevant fields in order to build the scaling theory. One 
possible strategy could example be to adapt the Dynamically Driven Renormalization Group (DDRG) procedure 
developed by Pietronero et al. [30]. We believe that exponents analogous to a, rj can be extracted for the Zhang model 
from the scaling of the expectation with respect to fi^ and of the spectral gap of W. 

For the generating function (59) of the marginal probability distribution Phin) the scaling form (60) reduces to : 

GL{t) = L<^-^^-^-^G{tL^-) (61) 

where G is a universal function. This form corresponds to the Kadanoff et al. finite-size scaling ansatz [23]. In 
particular, the cumulant CL{q) has a scaling factor : 

I ^ def log(CL(g)) 
(j{q) = hm 62 

L^oo log(_L) 

which is a linear function^^ of q, a{q) = /3ra((/ + 1 — t,,), for <? > r,,, — 1. 

Note that the measured exponents t„ G]1,2[ for the usual avalanche observables. Consequently, a{q) > for 
q > Tn — 1, and the moments of order g > t„ — 1 diverge in the limit L ^ oo. This implies a loss of analyticity of 
the limiting generating function. As discussed below and in [13] this effect is manifested by a Lee- Yang phenomenon. 
In particular, the properties of the critical zeros are directly related to the a{q)''s. In the case (61) the distribution of 
zeros is characterized by the exponents f3n,Tn, but relations linking the Lee- Yang zeros distribution to the cr(q')'s can 



'^^Recently, it has been argued that the finite-size form (61) might be violated in several models like the BTW model [31] or 
the Zhang's model [32]. Alternative scaling where a{q) is a non linear function of q have been proposed [31]. The conclusions in 
[31,32] were however essentially supported by numerical simulations and theoretical results are still missing. In particular the 
numerical bias induced by numerics were not discussed [13]. Such a scaling entails very singular properties for the asymptotic 
distribution function (resp. the asymptotic topological pressure, if it exists) and, consequently, for the asymptotic dynamical 
system. This opens interesting questions and perspective discussed in a forthcoming paper. 
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be derived for more general scaling [13] . 



Let us now discuss another extension of the potential (47). At a microscopic level, this potential characterizes the 
local contraction along a trajectory. Since the local contraction is related to the fractal structure of the support of 
the invariant measure, this remark suggests that it should possible to make a connexions between the multifractal 
spectrum of the measure jlL and the avalanches distribution. More precisely, the dynamical system (33) may be 
considered as a probabilistic graph iterated function system (see [33] for a review) with maps {Tu,}^^^j, with transition 
graph A and where the probability is the SRB measure [7]. Indeed, though the maps are only quasi contractions 
any finite composition of those maps along the graph is a contraction. In order to apply the theory developed for 
iterated function systems one needs additional conditions (open set condition, or strong separation condition) which 
are discussed in [7]. Were the Zhang's model a conformal iterated functions systems would the potential which allows 
to compute the multifractal spectrum be given by: 

cPl^iCo) = -q\og{N) + /31og(e)s(wi) (63) 

The multifractal spectrum is the Legendre transform of the function D{q) = where /? is such that [33]: 

^liMiU) = (64) 

Though the Zhang's model is not conformal since the contraction is not uniform in the phase space (see the Lyapunov 
spectrum Fig. 1), the potential (63) might be useful when considering the marginal energy distribution of one site. 
This is under current investigations and will be published elsewhere [34] . 

To compute the multifractal spectrum of fiL a more elaborate potential is required [37]. In particular, the cor- 
responding thermodynamic formalism is not additive, like in the appendix, but sub-additive. The framework for 
such a generalized thermodynamic formalism is described in [35] and in greater generality in [36]. Let ai{uj,k), i = 
1 . . . N,k = 1 ... 00 be the singular values of the map T^^o. . .oT^^ ordered such that 1 > ai{u),k) > ... > aN{u},k) > 
and remark that there is a finite k such that, whatever w, VZ > fc, 1 > ai{u}, I). Define : 

gi3{u),k)=ai{ui,k)a2{oj,k)aj-i{u),k){aj{Ld,k)f~^~^'^, if < (3 < N (65) 
5^ (d;,fc) = (dei(T„,o...oT<,J^) = e^n=i if (3 > N (66) 

where, in eq. (65), j is the integer such that j — 1 < /3 < j. Consider the the (sub-additive) potential: 

4,<i,k{^) = -MN) + log{g0{Q, k)) (67) 



and define the pressure by; 



= l™ Uog <PhA^) (68) 

K — ^OO A/ 



Then there exists a unique (3 = (3{q) such that the topological pressure vanishes [37]. The numbers D{q) = are 
the Reyni dimensions, and the multifractal spectrum is the Legendre transform of £>(</). We are currently investigating 
the relations between this potential and the corresponding multifractal properties to macroscopic transport properties 
[34]. In particular, we are trying to interpret the singular values in terms of avalanches properties. 
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III. SOME REMARKS ABOUT THE THERMODYNAMIC LIMIT. 



In this section wo want to discuss several effects observed when the size L tends to infinity ( "tliermodynamic 
limit"). These effects can be anticipated with the tools and c;oncepts developed in the previous sections, though they 
do not allow yet to handle them analytically. Consequently, most of the results presented are conjectured from the 
thermodynamic formalism description and are numerically checked. 



Rather than attempting to define the thermodynamic limit for the Gibbs state one can, as usual in equilibrium 
statistical mechanics of phase transitions, focus on the analyticity properties of the thermodynamic potential (free 
energy). More precisely, from eq. (56), the topological pressure is nothing but the generating function for the joint 
distribution of the avalanche observable n. This quantity certainly exists for finite L and has to be also defined in 
the thermodynamic limit if we admit that there exists a limit probability when L — > oo. However, if a critical state 
is indeed achieved, the finite size topological pressure should therefore develop singularities as L — > oo, /i = 0. 

The topological pressure (56) is a complicated object to handle, even numerically. However, one can argue that 
if the generating function (58) of the marginal distribution develop singularities in the thermodynamic limit, then 

develops singularities as well. In equilibrium statistical mechanics, a standard way to handle the singularities of 
the free energy and the scaling as L — > cx) is the study of the Lee- Yang zeros [14]. In many examples the partition 
function of the finite size systems is a polynomial in a variable z which typically depends on control parameters like the 
temperature or the external field. Since all coefficients are positive there is no zero on the positive real axis. However, 
in the thermodynamic limit, at the critical point, some zeros pinch the real axis at 2; = 1, leading to a singularity in 
the free energy. The finite-size scaling properties of the leading zeros and of the density of zeros near z = 1 determine 
the order of the transition [38] and also the critical exponents in the case of a second order phase transition [39]. In 
this paper, we show numerically that the same effect arises for the generating function of PL{s,h) for h = 0, while 
there is no Lee- Yang phenomenon for h> 0. This property is not specific to the observable s or to the Zhang's model. 
We showed indeed in [13] that this property arises as soon as a probability distribution Pi,(n), n = l.-.^E < 00 
converges to a power law n = 1 . . .00 when L — > 00. Furthermore, when 1 < r < 2 (which is the case for the 
usual avalanche observables and all the models), interesting anomalous finite size scaling effects are observed. 

Since ^£(/i) is finite for finite L, the generating function 



where z e C, is a polynomial of degree particular, since 2l{z) is an analytic function of z in the complex 

plane, all its moments exist. Since all PL(n,h) are positive Zl{z) has no zero on the positive real axis for finite L. 
Consequently the log-generating function log{ZL{z)) is well defined on R^. More precisely, Zl{z) > for 2; in a small 
neighborhood of where Gz,(t) = log{ZL{e*)) is an analytic function of z (resp. t). 



A. Lee- Yang phenomenon. 
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(69) 
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For finite L, 2l{z) has + 1 zeros in C, that are either real < 0, or complex conjugate. Denote them by 

ZL{k),k = O...C£ and order them such that < \zl{1) - 1| < ... < \zL{k) - 1| < ... < \zl{^2 + 1) - 1|- Write 

We report now the following numerical observations for the avalanche size distribution. A general theory and 
analytical results, extending to other models of SOC, are developed in [13] . 

• The zeros lie on a curve Cl in the complex plane which accumulate to the unit circle. This curve is however not 
a circle. In particular, it has a "cusp like shape" in the neighborhood of ^ = 1 (Fig. 5a). An analytic form is 
given in [13]. 



Lee-Yang zeroes in the z-complex plane. Lee Yang zeroes in tlie z complex plane 




FIG. 5. Lee Yang zeros in the z plane. Fig. 5 a. The curve Cl for L = 40. Fig. 5 b Local behavior near a = 1 of the 
Lee- Yang zeros for various L values in the z complex plane Ec = 2.2, e = 0.1, h = 0. 

• For h = infinitely many zeros accumulate on z = 1 (see Fig. 5 b and Fig. 6). Consequently, the pressure 
is not analytic any more in the thermodynamic limit, as expected. Note that this property is equivalent to the 
divergence of the moments mL{q) since it can be proved that the moments mL{q), for any q larger than some 
Qo diverge if and only if a Lee- Yang phenomenon occurs. This opens a way toward a scaling theory from the 
behavior of Lee- Yang zeros, in a way similar to what was done in equilibrium statistical mechanics [38^2]. 

• For h > (sub-critical regime) the zeros do not pinch z = 1. This is a clear indication that in this case there is 
no critical state in the thermodynamic limit (Fig. 6). 
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FIG. 6. Distance of the first Lee- Yang zeros zi^{l) to z = 1 versus L, Ec = 2.2, e = 0.1 in the critical case h — and 
subcritical case {h = 0.1, h = 0.05). 



• There exist characteristic features of the zeros that can be connected to the critical behavior. 

— For h ~ 0, the angle 0]^{k) (argument of Z[^{k)) scales like 6'i(fc) ~ (Fig- 7). It can be proved that, for 
a probability distribution converging to a power law and obeying a finite size scaling form s^'^g{^), this 
scaling is exact (with however a slight deviation for the 2 first zeros [13]). More generally, if the probability 
distribution obeys the scaling form (29) then one can prove that the angle 6'l(A:) is scales like: 

For n = s, the avalanche size, it follows therefore from (25, 30) that 6'i(fc) ^ -^^^ namely that 

0L{k) (71) 

[a + cL-^] " 

where a,c are respectively proportional to 2(1 — 'y)eX^ and (see eq. (25)). a, c depend slightly on 
L (via X^) but they converge rapidly to a constant as L — > cxd (see fig. 4a for the L dependence of 
X^). Therefore, studying the L dependence of 9L{k) gives a straightforward way to compute a (resp. 
(3 — ^). We performed a numerical study of the angles 0]^{k) for — 2.2, e = 0.1 in the conservative 
and non conservative case {h = 0.1), for i = 10 — 50. We obtained /3 ~ 2.59 ± 0.04 and cr = 0.772. This 
corresponds to a critical exponent t — 1.227, not so far away from the theoretical value r — 1.253 obtained 
by the renormalization group analysis [30]. Going to larger size would certainly improve the accuracy, 
though one has to be careful to pathological bias induced by the standard numerical procedure where the 
number of sample is fixed independently of the system size [13]. In Fig. 7 we plotted the angle ^^(fc) for 
fc = 4, 6, 8 and a non linear fit where the constants a, c in (71) where used as fitting parameters. We found 
a = 0.00439 ± 2.10-^ and c = 1.684 ± 3.10-3 
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FIG. 7. Argument of the Lee- Yang zeros 2^(4), 2^(6), (8) for Ec — 2.2, e — 0.1 in the critical case h = and in the 
subcritical case h = 0.1. (inset). In full color lines are drawn the interpolation curves y — 2Tvk{cx~'^)~ (conservative case) and 
y — 2Tvk{a + cx^^)~ (non conservative case) where a, c have been determined by non-linear fitting. 

— In many models of statistical mechanics the singular part of the free energy obeys a finite size scaling form 
f''{t,V) = yW[t{AV)^] (equivalent to (61)) [29,41,42], and the first zeros in the t plane (t = log{z)) 
do a characteristic angle with the real axis, independent of L, which can be related to the specific heat 
exponent a. Here the angle in the t plane violates the usual scaling and depends weakly on L (in fact the 
violation is logarithmic [13]). It can be proved that this effect is not an indication that the topological 
pressure does not obey finite size scaling but is simply due to the value of the exponent t > 1. 

B. Initial conditions sensitivity. 

The presence of a singularity set S for the dynamics induces an interesting phenomenon. If a trajectory approaches S 
sufficiently close at a time t say, a small perturbation in the energy configuration at time t can induce a response which 
is not proportional to the perturbation. This happens if the perturbed trajectory crosses the singularity set since the 
avalanche will be different in the initial configuration and in the perturbed configuration. This effect is obviously due 
to the presence of a threshold in the dynamics definition. More precisely, if in the trajectory of an energy configuration 
X a site i is such that its energy is arbitrary close to Ec at some time t, then obviously a small perturbation on this 
site can induce a completely different evolution. This phenomenon is particularly prominent if the measure of any 
?7-neighborhood of the singularity set is positive, for rj > 0. Indeed, in this case, with positive probability a trajectory 
will show non linear sensitivity to arbitrary weak perturbations for those times when it approaches the singularity set. 
In other words, in this case, taking two arbitrary close typical energy configurations and exciting the same sites along 
the whole trajectory of these two configurations, there will exist a finite time such that each configuration undergoes a 
different avalanche. This will result in an effective unpredictability of the evolution (weak initial condition sensitivity) 
and a failure of differentiability, or in more physical terms a failure of linear response, of the observables along typical 
trajectories. 

It is therefore instructive to compute numerically the probability that a site is close to Ec within a distance 
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lower than some t] > 0. Call Uri{S) = {X|3i, \Xi — Ec\ < rj}. Then fj.^ {Urj{S)} is exactly the probability that, 
during its evolution, a site is close to Ec within a distance lower than rj. In Fig. 8 we have drawn the variation of 
/iL as a function of Ec for L = 30 (Fig. 8a) and as a function of L (Fig. 8b), for rj — 0.001 and ?/ = 0.0001, 

Ec — 2.2, e — 0.1, h = Q. The statistics were done for 10 trajectories and 10^ time steps per trajectory. 
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FIG. 8. Probability hl {Ut^{S)} that, during its evolution, a site is close to Ec within a distance lower than some rj, as a 
function of Ec, for L = 30, e = 0.1, h^O (Fig. 8a) and as a function of L, for Ec = 2.2, e = 0.1, h = (Fig. 8b). 



In Fig. 9 we have plotted fi^ {^^jj('5)} as a function of r] for L = 30, Ec = 2.2, e = 0.1, h = 0,h = 0.1. Though 
we were very careful in computing this probability (increasing the number of samples when rj decreases) we believe 
that the left most points are biased. Ignoring this point we observed that fi^ {^^jj('5)} decay like 77", as 77 — > 0, where 
a — 0.98 ± 0.01 for h = and h = 0.1. This is supported by another well known result, already in the early paper 
of Zhang. The projected energy density on any site has a characteristic structure depicted Fig. 9 b. More precisely, 
in this figure we have plotted p(-E) fiL [X | 3i G A, Xi G [E — l, E + l]]. The remarkable peaks structure of this 
distribution has interested many people but nobody has been able, up to now, to give an analytic form for it. However, 
this is not the point that interests us. Rather, if we focus on a small neighborhood of Ec, we note that the p{E) is 
seemingly absolutely continuous on the left (explaining the exponent a ^ 1). 
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FIG. 9. a. Probability ^l{W^(5)} as a function of q for L = 30, = 2.2, e = 0.1, h = 0, ft = 0.1. Fig. 9 b. Probability 

p{E) =^ HL [X I 3i G A, Xie[E-t,E + l]], where d = 2, L = 30, = 2.2, e = 0.1, ft = 0, t = 10"^ Insert: zoom about 
E,. 

The conclusions are quite interesting. On one hand, this last result numerically validates the assumption (40) 
we made in section II C to infer the existence of a finite Markov partition. On the other hand, it shows that the 
probability to approach a small neighborhood of <S, though weak, is non zero. Consequently, the Zhang model displays 
this interesting form of initial condition sensitivity and linear response violation. It would be interesting to adapt 
this analysis to other models, in particular those relevant for earthquakes. Moreover, {Un{S)} decreases with Ec- 
Finally, fiL {UniS)} increases with L. This suggests that this phenomenon is more and more prominent as L grows. 
Note, that the number of singularity planes diverges as L ^ oo, for /i = 0. It might then be that the singularities 
become dense, leading to an extremely impredictible and singular system. 

IV. DISCUSSION AND CONCLUSION. 

In this paper we have discussed the application of thermodynamic formalism to the Zhang's model of Self-Organized 
Criticality. We have shown that under physically natural assumptions like ergodicity Gibbs measures can be defined to 
characterize various statistical properties of the finite size SOC state. This opens perspectives to build the equivalent 
of the statistical mechanics theory of critical phenomena for SOC models. It may open the way toward a general 
setting in which concepts like universality classes could be properly defined. Indeed, though the extrapolation of 
the thermodynamic formalism to the infinite lattice size limit needs further developments, this work suggests that 
the topological pressure can be used as an indicator of a phase transition. In particular we exhibited a Lee- Yang 
phenomenon for a quantity derived from the pressure and we noticed that several characteristic patterns emerges, 
which could allow a classification of the models in the future. 

We would like to discuss now some points that have not been developed in the paper. 

1. Extensions of the thermodynamic formalism. In this paper we have focused on the most common potentials 
which are directly related to dynamical properties and to the fractal structure of the support of the invariant 
measure. We also shown that they are related to the avalanche size. We now intend to construct more general 
potentials allowing on one hand to investigate the properties of other avalanche observables. On the other hand, 
we showed numerically in [8] that the Lyapunov spectrum exhibits a finite size scaling property with an universal 
exponent t\ related to the anomalous diffusion exponent. It would be worth to show this property analytically 
by producing a suitable potential. Finally, we would like to study the action of the dynamical renormalization 
group analysis developed in [30] on those potentials. 

2. Spectral gap vanishing. In the finite system, the exponential correlation decay along the time trajectory is given 
by the spectral gap between the largest eigenvalue of the Markov matrix W (which is 1) and the second eigenvalue. 
As mentioned above, this gap is positive whenever A is mixing. When L diverges, Ml diverges for /i = since 
A/j, > On the other hand, the notion of critical phenomena somehow involves non exponential correlation 
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decay or a divergence of the time correlation length which is the inverse of the spectral gap. Consequently, one 
expects that, for /i = the spectral gap vanishes in the thermodynamic limit. It might be useful to compute 
the spectral gap and particular its L dependence. This could be achieved from standard techniques on Markov 
chains [43,?], provided we have additional informations about the structure of W. Prom the analogy with 
critical phenomena we expect the gap to vanish like L'^ where r] plays the role of the exponent giving the spatial 
correlation decay in statistical mechanics. This could be a new exponent that might be related to the exponent 
T\ that we have exhibited in [8]. 

3. Explicit form of the energy density. In the section II C we used the result /xi, {Z//^(<S)} ~ that was only 
checked numerically. An analytic computation would be certainly be more satisfactory. It could be done if 
we could achieve the computation of the energy profile depicted fig. 9b This would also allows to elucidate 
its particular peaks shape which is still an open problem in SOC. This could be achieved from the transport 
equation developed in the section IB. 

These points are under current investigations. 
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APPENDIX 

In this appendix, devoted to non-specialists, we give a brief summary of the thermodynamic formalism used in 
section IID to construct Gibbs measures. Useful references are [12,11,16,48] 

Assume that we have a set of symbols u G I and a transition matrix A. Call X (resp. the set of bi-infinite 
(resp. right infinite) legal sequences. Write [u)]n for a n-cylinder (this is the subset of X where the sequences have the 
same n first digits as u). Denote by the shift on X. In the sequel we restrict to the set of right infinite sequence 
X~^ as usual in the frame of thermodynamic formalism [16,48]. There exists a formal analogy between a sequence 
of X~^ and a (right-infinite) one dimensional chain of Potts-like spins taking values in T. The transition matrix A 
acts then as a hard-core like potential in the sense that, if the spin at the t th place in the chain has a value Wt, 
the next spin (at the place t + 1) can only take values in the subset of T such that A^^.^t^i = 1. Other transitions 
are forbidden. This formalism allows in particular to study a class of dynamically relevant invariant measures of the 
dynamical system as formal analogous to Gibbs measure in a chain of spins. 

Let F{X~^) be the space of Holder continuous functions, for a metric d0{x,y) = 9^ on where A'' is the largest 
non-negative integer such that Xi = yi,i < N and < ^ < 1 [48,16]. We call a potential an element of F{X~^). A 
potential has in particular the following property: 

varn{(l)) < C6I", n > (72) 
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for some C > 0, some 9 s]0, 1[, where varn{4') = sup{|0(a;) — (j){y)\, xt =yi,i< n}. Note that eq. (72) is the 
equivalent of the exponential decrease of the interaction with the distance, insuring the existence of a thermodynamic 
limit in statistical mechanics [24,25]. A finite range potential of order r is such that (p{x) = ^{y) if Xi = yi for 
< i < r, namely the values that (j) takes depend only on the first symbols (resp. the r first spins). Any infinite range 
potential in can be uniformly approximated by a sequence of r-range potentials. 

Set S't?!>(w) = Yld=i ^i^A^) where ^ is a potential. Define the finite partition function by 

M<P) = Yl expSTcl>{u>) (73) 

and the pressure (resp. the free energy per spin) by : 

n'P) = lim ^logZri^) (74) 

I — *oo 1 

It can be proved that this limit exists provided (j) decays sufficiently fast (see eq.(72)), namely there is no "phase 
transition" provided the potential belong to F{X~^). 

A Gibbs measure is an invariant measure where the potential gives an exponential weight to the cylinders. More 
precisely, /x^ is a Gibbs measure if there exists a A such that for all T > and u € Xr 

This essentially means that ji^ is exponential with a weight given by the sum of the values that (j) takes on the orbit 

of ui. Since ZT{(j)) ~ e^^^'f') the mcasiirc of a spin chain of length T is ~ — ^^^Zrit) ~ ^'^'^ formal analogy with 

statistical mechanics is straightforward. 

In this setting one also associates to each potential (j) the Ruelle operator C^, which is a formal extension of the 
Kramers- Wannier transfer matrix for a spin chain [24]. An extension of the Perron- Frobenius theorem for matrices, 
due to Ruelle, shows that when A is irreducible and aperiodic admits a unique maximal positive eigenvalue 
which is equal to the pressure J-'{4>). The corresponding left eigenvector is the Gibbs measure /i^. The spectrum of 
C4, provides informations about the (strong) mixing properties of /i^. In particular, the spectral gap between the 
largest eigenvalue (J-'((j))) and the remaining part of the spectrum determines the dominant exponential decay rate of 
correlations fmictions (resp. decay rate to equilibrium). 

/i0 satisfies also a variational principle analogous to the free energy minimization in statistical mechanics. Namely, 
call h{p,) the entropy of the invariant measure fi then the quantity h(j.i) + J cpdfi admits a unique maximum for /i = /i^, 
equal to the pressure J-{(t>). The maximizing measure is naturally called an equilibrium state. Each equilibrium state 
for a potential G F{X^) is a Gibbs state. There exists only one maximum (resp. one equilibrium state related 
to the observable cp) when A is mixing. This situation corresponds to the absence of phase transition in statistical 
mechanics. 



The pressure, beyond the variational principle, shares others characteristic with the free energy: it is convex, non 
decreasing, sub-additive {J-{(pi + ^2) < •^('^1) + ^{4"2)) and this is a generating function for the expectations with 
respect to /x^. More precisely, the following can be proved [48,16]. Let ry G Fg, and set 'P(t) =' J-{4) + ti]), where 
i e R, then : 
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V'{0)= f^j^i<l> + tv) 



= f iid^i^''^' E[ri\^ (76) 

t=o J 



where £^[]^ is the expectation with respect to /z<^. In the same way : 

T 



= = ^Y.^ W)m]4, - E Ml (77) 

where E ['nit)r]{0)]^ — E [ry]^ is the correlation function of the function rj, at time t, the average being performed with 
respect to /x^. This can be generaUzed to correlation functions between different observables and to higher order 
[45]. The coefficient <J,f,{r]) characterizes the average ffuctuations of r] along trajectories weighted by the measure 
Provided cr^(?7) > the central limit theorem holds, namely the ffuctuations are Gaussian (more precisely 



limT->cx) Prob 



= ^f{y) where M is the characteristic function of the Gaussian distribution). 



Prom the relation (77) one can extract Green-Kubo transport coefficients from microscopic quantities [45,44]. 

According to the choice of the potential one is able to extract different informations about the statistical properties 
of the dynamics. The situation is somehow analogous to statistical mechanics where the choice of the thermodynamic 
potential corresponds to a different choice of ensemble. However one has a priori an infinite number of choice for the 
potential (resp. measure) but only a few of them are physically relevant. 
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